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1) Advection equation and method of characteristics. 
1.1 Advection equation. 

• We consider a given real number a > and we wish to solve the 

so-called advection equation of unknown function u{x, t) : 

/. . . N du du ^ „ ^ 

1.1.1) — + a— =0, t>0, xeM. 

at ox 



We first look to the homogeneity coherence of the different terms of equation 
(1.1.1). On one hand, the ratio ^ is homogeneous to the dimension [u] of 
function •) divided by the dimension [t] of the time and we have : ^ ^ 

On the other hand the expression a is homogeneous to the dimension 
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[a] of scalar a multiplied by the ratio and we have a ^ ~ [a] 
equation (1.1.1), the two previous terms ^ and have the same dimension 

and we deduce from the previous formulae the equality : ftj ^ • Then we 
have established that the constant a is homogeneous to a celerity : 

(1.1.2) N - ^- 

• The Cauchy problem for the model equation (1.1.1) is composed by the 
equation (1.1.1) itself and the following initial condition : 

(1.1.3) u{x, 0) = uq{x) , a: € IR, 

where IR 9 a: i — )■ uq{x) G 1R is some given function. We observe that the 
solution of equation (1.1.1) is constant along the characteristic (straight) lines 
that satisfy the differential equation 

(1.1.4) - = a. 

Proposition 1.1. The solution is constant along the characteristic lines. 
Let < A < t be some given parameter and ii(», •) a solution of equation 
(1.1.1). Then function w(«, •) is constant along the characteristic lines, i.e. 

(1.1.5) u{x — aX, t — X) = u{x, t) , \/ X, t, X. 

• The proof of Proposition 1.1 is obtained as follows. We consider 
a fixed point {x, t) in space-time IR x [0, +oo[ and the auxiliary function 
[0, t] 3 X I — > v{X) = u{x — aX, t — X) . We have, due to the usual chain rule 
for derivation of operators 
dv 



dX 



du du 

{—a}— — {x — aX,t — X) = if function n(«, •) is solution 

of the advection equation (1.1.1). Then ^;(A) does not depend on variable A 
and we have in particular v{X) = v{0) , which exactly expresses the relation 



3 



FRANgois Dubois 



(1.1.5). We have in particular for X = t : u{x, t) = u{x — at, 0) = UQ{x — at) 
as illustrated on Figure 1.1. □ 

t 





1 


J u(x,t) 






^ 


x—Xa. ^ 





UgCx-at) 

Figure 1.1. The solution u{x,t) of the advection equation 
is constant along the characteristic lines. 




1.2 Initial-boundary value problems for the advection equation. 

• The second step is concerned by the so-called initial-boundary value 

problem considered for x > and t > with some given initial condition 
Uo{x) for t = and a boundary condition Vo{t) for x = : 
du du 

(1.2.1) -— -|- a 7— = 0, t>0, X > 0, (equation) 
at ox 

(1.2.2) u{x, 0) = Uq{x) , x>0, (initial condition) 
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(1.2.3) u{0^ t) = vo{t) , t>0, (boundary condition). 

Proposition 1.2. Advection in the quadrant a: > and t>0. 

We suppose that a > 0. Then the solution of the advection equation (1.2.1) 
with the initial condition (1.2.2) and the boundary condition (1.2.3) is given by 
the relations 

(1.2.4) u{x, t) = uo{x — at) , x — at > 

(1.2.5) u{x, t) = vo(^t — , X — at < . 

The initial condition wo(«) is advected towards space-time point [x^ t) when 
X — at > and the boundary condition vo{») is activated for x — at < . 

• Proof of Proposition 1.2. 

In order to solve the problem (1.2.1)-(1.2.3), we use the method of characteristics. 
We fix a point [x, t) of space-time domain that satisfies a; > 0, t > and we 
go upstream in time with the help of the characteristic line that goes through 
this point (see Figure 1.2) : 

(1.2.6) x{X)=x-aX, t{X)=t-X. 

• First case : x — at>0. When we take the particular value X = t in the 

previous relation (1.2.6), the particular point y = x{t) = x — at on the axis of 
abscissa is strictly positive then the initial condition uo{y) is well defined. The 
solution w(», •) is constant on the characteristic line (see Proposition 1.1) that 
contains this particular point. Then relation (1.2.4) is established. 

• Second case I X at < . We consider the particular value A = - 
inside the expression (1.2.6). Then the corresponding foot of the characteristic 
belongs to the time axis : 6 = t — X = t— ^ and ^ > due to the inequalities 
X < at and a > . The solution is constant along the characteristic line going 
through this point and the relation (1.2.5) is established. □ 

• In the particular case where datum uo{x) is identically equal to zero, 
i.e. 

(1.2.7) uo{x) = , x>0, 

and if the boundary condition Vo{t) is sinusoidal for time positive to fix the 

ideas, 

(1.2.8) vo{t) = sin(wt) , t > 0, 

the solution of the advection equation in the domain a: > , t > via the 
relations (1.2.4) and (1.2.5) can be considered with the two following view points. 

(i) We take a snap shot of the solution m(», •) at a fixed time T > 0. 

We consider the partial function [0, -|-oo[ 3 x i — > u{x, T) e IR and taking into 
account the relations (1.2.4), (1.2.5), (1.2.7) and (1.2.8), we have 

(1.2.9) uix,T) = ^ ^(^-?) 



X < aT 



a 

, X > aT 
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and this function is illustrated on Figure 1.3. 

(ii) We fix a particular position X in space and we look, as time is 

increasing, to the solution •) at this particular point. We show on Figure 
1.4 the function [0, +oo[9 t i — > u{X, t) e IR and taking into account the 
relations (1.2.4), (1.2.5), (1.2.7) and (1.2.8), we have 

, t < 

(1.2.10) u{x, T) 



sm 



a 



u(x, T) 




Figure 1.3. Snap shot of the solution of the advection equation 

at time t = T. 

u(X, t) 




Figure 1.4. Evolution of the solution at the particular point x = X. 

1.3 Inflow and outflow for the advection equation. 

• We still suppose that celerity a is positive and we consider the resolution 
of the advection (1.2.1) in the space-time domain 
(1.3.1) < a; < L, t > 0. 

The relations (1.2.4) and (1.2.5) can still be applied because the proof of Propo- 
sition 1.2 remains unchanged in this particular case. As a consequence of the 
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previous property, we remark that no boundary condition is necessary at the 
particular position x = L for solving the advection problem in the space-time 
domain defined in relations (1.3.1). The initial condition (1.2.2) has simply to 
be restricted in domain ]0, I/[ : 

(1.3.2) u{x,Qi) = uq{x) , 0<x<L, 

and the boundary condition (1.2.3) at x = remains unchanged : 

(1.3.3) w(0, t) = vo{t) , t > 0. 



the boundary x=0 
is an input 
when a > 







the boundary x=L 
is an output 
when a > 



Figure 1.5. Initial-boundary value problem for the advection equation 
with a > in the domain < a: < L and t > 0. 



output 




u(L, e) = \\^(e) 



u(y, 0) = uj;y) 

Figure 1.6. Initial-boundary value problem for the advection equation 
with a < in the domain < a; < L and t > 0. 



• The difference between point a; = and point x = L for the resolution 
of the advection equation in space-time domain (1.3.1) is due to the fact that we 
choose an orientation of the characteristic lines x — at = constant associated 
to an increase for the time direction. With this choice of time direction, the 
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characteristic lines enter inside the space-time domain (1.3.1) at x = and 
they go outside at x = L . The boundary condition (1.3.3) is given at the input of 
the domain (see Figure 1.5) and at x = L , there is a free output from space time 
domain (1.3.1), without necessity to specify any numerical boundary condition. 

• If we change the sign of celerity a, i.e. if we suppose now 

(1.3.4) a<0, 

the above analysis remains unchanged, but the algebraic relations (1.2.4) and 

(1.2.5) have to be modified (see Figure 1.6). We still start from relation (1.1.5) 
that expresses that the solution of the advection equation (1.1.1) is constant 
along the characteristics lines. The foot of the characteristic line that contains 
the particular point {x, t) in space-time is either the point {y = x — at, 0) if 
X — at < L , either the point (L, 9 = t — ^{x — L)) if x — at > L . In the first 
case, we have y > and ^ < then the initial condition (1.3.2) is advected 
inside the domain (1.3.1) and we have : 

(1.3.5) u{x, t) = uo{x — at) , x — at < L . 

• On the contrary, if x — at > L, we have y > L and ^ > then the 
boundary condition at x = L that takes now the expression 

(1.3.6) u{L, t) = WL{t) , t > 0, 

is advected inside the domain of study and we have : 

t -\ ~~)' X — at > L . 

We have established the following 

Proposition 1.3. Advection in the domain 0<a;<Z/, a<0. 
Under the hypothesis (1.3.4), the resolution of the advection equation (1.2.1) 
in the space-time domain (1.3.1) conducts to a well posed problem when we 
introduce the initial condition (1.3.2) on the interval ]0, L[ and the boundary 
condition (1.3.6) at the input region located at x = L, without any boundary 
condition at the output located at x = 0. The solution of Problem (1.2.1), 
(1.3.2) and (1.3.6) is given by the relations (1.3.5) and (1.3.7). 



2) Finite volumes for linear hyperbolic systems. 
2.1 Linear advection. 

• We still study the advection equation parameterized by some celerity 

a > : 

(2.1.1) — - + ^{aW) =0, t>0, xeM, 

and we search a discrete version of this mathematical model. For doing this, we 
introduce a space step /S.x > and a space grid composed by points Xj whose 
coordinates are multiples of this space step Ax, id est 
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(2.1.2) 



j Ax 



For a finite domain, ]0, Z/[ to fix tlie ideas, tlie above grid is limited to integer 
values j such that 

(2.1.3) < , < J = A 

and the vertices {xj)o<j<j are usually used in the context of the finite differ- 
ence method. The intervals Kj_^i/2 = ]xj, Xj-^-i[ between two vertices can be 
considered as finite elements (or finite volumes in our study) and they cover the 
entire domain ]0, L[: 

(2.1.4) [0, L] = y [x, , x,+,] , 

0<j<J-l 

as proposed in the general context of meshes (see e.g. Ciarlet [Ci78]). We 
introduce also a time step At > and the discrete time values at integer 
multiples of the above quantum : 

(2.1.5) e = nAt, n e IN. 

We consider now a space-time volume Vj^_^^^(^ obtained by cartesian product 
of the two intervals ]xj , Xj-\-i[ and jt"', t'^~^-^[ (see Figure 2.1) : 



(2.1.6) 



V 



n+l/2 



]x,,x,+,[x]r,r+'[ 



i+i/2 



• The finite volume scheme consists simply in integrating the advection 

equation (2.1.1) inside the space-time domain V'^^^j^ introduced previously : 

dW d 



(2.1.7) 



n+l/2 



i+1/2 



dt 



+ 



dx 



{aW) 



dx dt 



, < j < J 



n > 



n+1 




n+1 



>1 ^^j-1/2 -j+1/2 -j+1 

Figure 2.1. Space-time grid for the finite volume method. 



Proposition 2.1. Finite volume scheme. 

Let IR X [0, -|-oo[9 {x, t) i — > W{x, t) e IR be a solution of the advection 
equation (2.1.1). We introduce the space mean value ^^^1/2 of ^^lis solution 
W{», •) in the cell '■ 
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(2.1.8) 



K 



J + l/2 




and the time mean value of the so-called flux aW{»^ •) at the space 

position Xj and between discrete times and t^^^ : 



(2.1.9) /, 



n+l/2 



1 

At 

Then we have the following constitutive relation of finite volumes schemes 



= a7 L {aW){xj,t)dt. 



(2.1.10) 



n+l 



-w. 



) + 



/ .n+l/2 _ nn+l/2\ 







This numerical modelling characterizes the so-called finite volume method which 
has been developed thanks to the work of S. Godunov [Go59], Godunov et al 
[GZIKP79], Patankar [Pa80], Harten, Lax and Van Leer [HLV83] or Faille, Gal- 
louet and Herbin [FGH91] among others. 

• The proof of Proposition 2.1 consists in a precise evaluation of the 

left hand side of equality (2.1.7). We use Fubini rule for the computation of 
double integrals and we begin by integrating in time for the ^ term : 



dW 



y^^^l'^ dt 



dx dt = 



J + l/2 



nX^l 



^Xj+i 



W{x, r+i)-T^(a:, r) 



dW 
~dt 



{x, t) dt 



dx 



dx 



Ax 



due to the definition (2.1.8). We proceed in an analogous way with the ^ 
term and begin now the Fubini procedure by integrating in space ; we have 



.n+l/2 
' J + l/2 



^Xj+l Q 

— {aW)ix, t) dx 
ox ^ ^ 



dt 



/ [{aW){xj+i,t)-{aW){xj,t) 



dt 



At 



pn+l/2 



n+l/2 



according to the definition (2.1.9). We add the two previous results, use identity 
(2.1.7) and divide by At Ax. We obtain exactly the relation (2.1.10). □ 

• The relation (2.1.10) is a very general form for the evolution of the mean 

values between two time steps. The increment {W^^^^^ ~ ^J+1/2) 

is, after correction by a multiplicative factor, equilibrated by the fiux difference 

(/ 



n+l/2 



i+1 ^ 



n+l/2 



) . The idea of a finite volume scheme is to consider now that 
the algebraic object is nomore the mean value of the exact solution 

but an approximation of this mean value. Then the relation (2.1.10) proposes 
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a numerical scheme for the discrete evolution of the approximated mean values 
7 j = 0, • • • , J— 1. Nevertheless, the numerical scheme is not entirely 
defined by the relation (2.1.10). Starting from mean values at the initial time 
step, i.e. 

1 f^^^ 

(2.1.11) = ^ Wo{x)dx, j = 0,---,J-l, 

we are able to increment the time step with relation (2.1.10) only if all the fluxes 
j.n+i/2 ^ j — Q . . . J j^ave been a priori first determined as a functional of the 
previous values. In a very general way, we say that the finite volume scheme 
(2.1.10) is an explicit scheme if each flux [q ^ given function of 

the mean values {W^_^-^^^) at the preceding time step number n : 

(2.1.12) /;+^/' = ^,({l^,V/2, fc = 0,---, J-1}), ^ j = 0, •••,J-1. 
The function vE'j is called the local numerical flux function at point xj and, 
joined with the evolution equation (2.1.10), its choice determines the numerical 
scheme. 

• A natural hypothesis claims that we have translation invariance for 

the evaluation of the flux if we move the discrete data in the same way ; in other 
words, the numerical flux function only depends on the p first neighbors 
of the interface Xj . Then the explicit numerical flux is a given function $ of 
the p first neighbors and we have : 

(2.1.13) /;+^/^ = ^(VFjV,/,^,, ■ ■ • , W^_,/„ WjVv^, • • ■ , • 
A very important particular case is one of a two-point scheme for the evaluation 
of the numerical flux. We have in this particular case : 

(2.1.14) /;+^/^ = $(i^;_i/2,^;+i/2)- 

With this particular choice, the numerical scheme for incrementing in time of 
the mean values takes the form : 

ii(w'Av2-w'iVi/2) + 

+ ^(*(W?+l/2. 1^^3/2) -*(WS"-l/2.W;+l/2)) = 0. 

It is also a three-point finite difference scheme. The finite volume scheme (2.1.10) 
(2.1.13) is said to be consistent with the advection equation (2.1.1) when the 
numerical flux function $ satisfles the condition 
(2.1.16) ^(W, ■■•,W) = aW , VVFelR. 

• The crucial question is how to choose a numerical finite volume scheme. 
The simplest choice consists in a two point explicit scheme such that the fi- 
nite difference scheme is identical to the upstream-centered scheme (see e.g. 
Richtmyer-Morton [RM67]). It takes the following expressions : 
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(2.1.17) ^ - W?+i/2) + «(wjVi/2 - W7-1/2) = , a > 

(2.1.18) i-(w;«,-w;Vi/2) + «(w';+3/2 - =0, a<0. 

The corresponding flux function is called the first order upstream- centered 
fiux, is simply given by the following relations : 

(2.1.19) ^W,W.) = [l^l'^ «>0 

When this flux function acts at a given point Xj of the mesh, we have : 

(2.1.20) f^" = #(W'J^,/„ W'Ai/2) = { awfZ, a < 0. 

If a > 0, the exact solution of the advection equation propagates the information 
from the left to the right ; the flux at the interface Xj is issued from the cell 
at the left of the interface and this cell at the number j — 1/2. If a < 0, the 
propagation of the information with the advection equation is from right to left ; 
the interface flux at the abscissa Xj is due to the control volume on the right, 
i.e. with number j + 1/2 as depicted on Figure 2.2. 



W 



j-1/2 



a>0 



a<0 



W 



j+1/2 



X. 

J 



Figure 2.2. Upwinding of the information for the advection equation. 





W, 



1/2 




x„=0 



j-l/2 



l-= 


=H 


Pn 


pH 


— 1— ^ 




1 
1 


\ ■ 




X. 

J 


- Vi 


1 
1 




Figure 2.3. Notations for the one-dimensional finite volume method. 
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• Recall that practical use of the upwind finite volume scheme like (2.1.17) 
when a>0 or (2.1.18) if a < is restricted to the usual Courant-Friedrichs- 
Lewy stability condition : 

(2.1.21) < 1 

as developed e.g. in the book of Richtmyer and Morton [RM67]. 
2.2 Numerical flux boundary conditions 

• In this section, we focus on the problem of the numerical boundary 
conditions. Recall that we study the advection equation in the space domain 
[0, L\ : 

(2.2.1) < a: < L 

and = € IN control cells (or finite elements) have been used to define a 
mesh : 

(2.2.2) JAar = L. 

Note that the cell is exactly the interval Xj\ and it is centered at 

point Xj_\i2 as shown on Figure 2.3. 

• At time step n At , the discrete field is entirely known and is composed 
of all the values W^_-^^2 j = " ' ^ J ■ With a flux function $(•, •) as 

proposed at relation (2.1.14), we observe that the two boundary fluxes f^^^^"^ 
and y^+^Z^ are not a priori deflned because states W^^^^ ^j+1/2 does not 
exist. The situation is more complex with numerical fluxes that use four points 
or more as proposed in (2.1.13) and will not be detailed in this section. Even 
if the formula giving the numerical flux at the boundaries has to be speciflcally 
studied, the finite volume scheme remains defined by the relation (2.1.10) and 
we have for the two cells encountering the boundary : 

(2.2.3) ^ (Wf+i - Wl),) + £ i/r"' - C") = , 

(2.2.4) (w^;« , - wu,.) + ^ ur''' - /ry') = o . 

• The question is now to adapt the relation (1.2.14) in order to determine 

the two boundary fluxes /(^^^^'^ at the left of the domain and at the 

right. For the advection equation with celerity a > 0, we have observed in the 
first section that some boundary condition vo{t) has to be assigned at a: = 
and it is not the case for x — L. It is therefore natural to take into account this 
information at the input of the domain and to set : 

(2.2.5) /o^+^/^ = ^ [av,{t)dt 
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or simply 

(2.2.6) /o"+'/' = avo{{n+^)At), a>0, 

if function t i — y vo{t) has a slow time variation at the scale defined by the 
time step. At the output x = L, no numerical datum has to be assigned to set 
correctly the continuous mathematical problem. We must maintain this property 
if we wish the numerical method to follow the mathematical physics as efficiently 
as possible. A simple boundary flux is associated with the previous numerical 
upwind scheme. For x = xj = L and a > 0, we observe that the upwind 
scheme (2.1.20) is simply written as : 

(2.2.7) = , a>0, 

and this relation (2.2.7) defines a first order extrapolated boundary flux. 

• The roles are reversed when a < 0. The abscissa x — corresponds to 
an output for the advection equation and the right boundary x = L is an input 
where a time field t i — > WL{t) is given. In the first case, the upwind scheme 
(2.1.20) can be applied without modification : 

(2.2.8) /o"+'/' = aW^/^, a<0, 

and it corresponds to a first order extrapolation of the internal data 

j = 1, • • • , J } at the boundary at time step nAt. For x = L, the boundary 

flux uses the given information between the two time steps : 

(2.2.9) /;+'/^ = awL({.n+\)At) , a<0. 

Proposition 2.2. 

Flux boundary conditions for the advection equation. 

When we approach the advection equation (2.1.1) with the finite volume method, 
the numerical boundary conditions induces a choice for the two boundary fluxes 
jn+i/2 jn+i/2 -^j^gj^ a > 0, the boundary Condition VQ{t) at the input 

can be introduced into the boundary with the relation (2.2.6) and the free output 
at the right can be treated with an extrapolation of the type (2.2.7). When 
a < 0, the free output at the left of the domain can be taken into account with 
the help of relation (2.2.8) whereas the input condition WL{t) at the right can 
be introduced thanks to relation (2.2.9). 

2.3 A model system with two equations 

• Let a > and 6 > be two positive real number. We study in this 
section a model problem that is composed by the juxtaposition of an advection 
equation with celerity a and an advection with celerity —h. We explicit the 
associated algebra : 

du du 

(2.3.1) — + a— = 0, t>0, xeJR, 
at ox 
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dv , dv ^ 

(2.3.2) 6— = 0, t>0, xeM. 

at ox 

We associate the two equations (2.3.1) and (2.3.2) and consider a unique problem 
with a vector field as unknown. We set : 

(2.3.3) V = (l) 

and the set of equations (2.3.1)-(2.3.2) can naturally be written as a system : 

S - (o -I) % - «■ 

By introducing the flux function -F(v?) according to the relation 
(2.3.5) F(^) = 

the system (2.3.4) takes the general conservative form : 

I + ^(^(^)) = «• 

• The approximation of system (2.3.6) with a grid parameterized by a 
space step Aa: and a time step At is conducted exactly as in the case of the 
advection equation. The following property is a straightforward generalization 
of Proposition 2.1. We left the proof to the reader. 

Proposition 2.3. Finite volume scheme. 

Let IR X [0, +oo[ 9 (a:, t) i — > (p{x^ t) e IR x IR be a solution of the linear 
conservation law (2.3.6). We define the space mean value ^^^^ solution 

(/?(•, •) in the cell : 

(2.3.7) ^-+1/2 = / ^{x, n dx 

and the time mean value of the fiux function introduced in (2.3.5) at 

the space position Xj between discrete times and t^^^ : 

(2.3.8) f-+^l^ = l^j^^ F{v(xi,t))At. 

We have the following relation that characterizes the finite volumes schemes : 

• We have now to propose a precise numerical flux function analogous 
to the relation (2.1.12) to transform the conservation property (2.3.9) into a 
flnite volume numerical scheme able to propagate the discrete values ¥'j_|_i/2 
up to the discrete time For internal interfaces xj, j = 1, • • • , J— 1 , it is 
natural to apply the upwinding scheme (2.1.20) with a left upwinding for the flrst 
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equation and a right upwinding for the equation (2.3.2). Figure 2.4 illustrates 
the associated algebra : 



(2.3.10) /; 



n+1/2 



u 



n 



a>0 



j-1/2 



V. 



j-1/2 



u 



j+1/2 



-b<0 



n 

^j+1/2 



X 



Figure 2.4. Interface upwind numerical flux 
for a model problem with two equations. 



u„(t) 



® 




Figure 2.5. 



x=0 x=L 
Boundary conditions for a model problem with two equations. 



• At the left boundary a: = , we have an input for the variable u and 

we suppose given the associated datum [0, +oo[9 t \ — > uo{t) e IR : 

(2.3.11) w(0, t) = uo{t) , t>0 

whereas it is an output for the v variable. By association of relations (2.2.6) 
and (2.2.8), we obtain 

(2.3.12) fr'^ = {"""^j^^*^) ■ 

At the other boundary of the interval ]0, I/[ , we have an output for the first vari- 
able u and an input for the second one, and an associated boundary condition 
[0, +00 [9 1 1 — y VL{t) is supposed to have been given : 

(2.3.13) v{L, t) = VL{t) , t > 
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as illustrated on Figure 2.5. The numerical flux at the right is evaluated by 
association of the relations (2.2.7) and (2.2.9) : 



2.4 Unidimensional linear acoustics 

• We consider a gas in a pipe of uniform section at normal conditions 
of temperature and pressure. The reference density is denoted by po and the 
reference pressure is named po- The sound celerity cq of this gas satisfies the 
relation 

(2.4.1) CO = .1^ 

V Po 

with 7 = 1.4 as proved e.g. in the book of Landau and Lifchitz [LL54]. A 
sound wave is a small perturbation of this reference state. The differences of 
density, pressure and velocity fields are denoted respectively by p, p and u. The 
hypothesis of a small perturbation implies that the entropy of the reference state 
is maintained for all the time evolution and in consequence, it is easy to establish 
the following relation between the perturbations of density and pressure : 

(2.4.2) p = clp. 

• The conservation of mass leads to a first order linear conservation law : 
, dp du 

and the conservation of momentum links the time evolution of velocity with the 
spatial gradient of pressure : 
. .^ du dp 

We introduce the vector 1^ = ^ of unknowns. Then the equations (2.4.3) 

and (2.4.4) can be written as a linear hyperbolic system of conservation laws : 
dW ^dW 

with 

(2.4.6) A=(l "f). 

• When we consider the eigenvalues and eigenvectors of matrix ^4, it is 
natural to introduce the characteristic variables defined respectively by 

(2.4.7) = p + pqCqU 

(2.4.8) (f- = p - pqcqu 

and the quantity po cq is named the acoustic impedance. We have from the 
relations (2.4.3) and (2.4.4) : 
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diOj^ d(p-L /dp du\ / dp 2 

-W + '^^ = {di + P'"'ai) + + 9iJ 

2 fdp du\ / du dp\ 

dip- dip- fdp du\ fdp du\ 

-W - = Kdi- di) - "° iai - """^ d^) 

2 fdp du\ f du dp\ 

and we recover a system of the type (2.3.4) studied previously : 

• A typically physical problem is the following : a given acoustic pressure 

wave [0, +00 [ 3 t 1 — > n(t) > is injected at the left a: = of the pipe and 
the waves go away freely at the right boundary x = L . At t = 0, the velocity 
and pressure of the fluid are given : 

(2.4.10) u{x, 0) = uq{x) , < a; < L 

(2.4.11) p{x,Qi) = pq{x) , 0<x<L. 

From a mathematical viewpoint, the boundary conditions have to respect the 
dynamics of this system of acoustic equations written in diagonal form (2.4.9) : 
the variable must be given at a: = and the variable (p- at the abscissa 
x — L. From (2.4.7) and (2.4.8), we determine the pressure as a function of the 
two characteristics variables (^+ and (p- : 

(2.4.12) p = 

and if the pressure is imposed at a: = 0, the relation (2.4.12) can be written 
under the form : 

(2.4.13) ^+{0,t) = -(^_(0, t) + 2n(t) , x = 0, t>0, 

that makes in evidence a reflection operator : the input variable is a given 
affine function of the output variable (p- . At the other boundary x — L , the 
notion of free output expresses that the waves that go outside of the domain 
of study have no reflection at the boundary. When x = L ^ the characteristic 
variable is going outside and there is no boundary condition for this variable. 
We have to express also that this wave has no influence on the characteristic (p- 
that wish to go inside the domain ]0, L[. In other terms, the input value p- is 
independent of the variable and also of time. We have in consequence 

(2.4.14) §i^-iL:t) = 0. 
We have established 
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Proposition 2.4. Boundary conditions for acoustic problem. 

The mathematical boundary conditions associated with the datum of a given 
acoustic pressure wave [0, +00 [ 3 t 1 — > n(t) > at the left of the domain 
]0, L[ admits the expression (2.4.13) and a condition of free output of the waves 
at the right boundary x = L can be expressed by the relation (2.4.14). 



P = n(t) 




1 



cp (L , t) = 9 (L , 0) 



Figure 2.6. 



Solution of the acoustic equations in one space dimension 
for a model problem with two equations 



• The above acoustic problem associated with the first order partial dif- 
ferential equations (2.4.3) (2.4.4), the initial conditions (2.4.10) (2.4.11) and the 
boundary conditions (2.4.13) (2.4.14) is illustrated on Figure 2.6. The initial 
conditions are active in the beginning of the evolution in time (t < ^) and 
have a trace for higher times due to the boundary conditon (2.4.13), that gives, 
due to (2.4.8) and (2.4.13) : 

(2.4.15) (p-{x,t) = p{x, t) - poCou{x, t) = po{L) - po Co uo{L) , t>—. 

Co 

On the other hand, the inflow boundary condition (2.4.12) and the second row 
of matrix equation (2.4.9) implies : 

r (p+{x, t) = p{x, t) + PqCqU{x, t) = 

^ ^ 1 = 2n(t--) - (^_(o,t--), t>-. 

K Co Co Co 

We deduce from the relations (2.4.15) (2.4.16) joined with the definitions (2.4.7) 
and (2.4.8) : 

(2.4.17) p{x,t) = U{t-fJ, 0<x<L, t>^^ 

(2.4.18) u{x,t) = uo{L) + ^^{li{t-^)-po{L)), 0<x<L, t>^. 

• We turn now to the numerical finite volume scheme. We have to de- 
termine the internal fluxes f"^+-^/^ ^ ji = l,---,J— 1 and the boundary fluxes 
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(2.4.21) /, 



jn+i/2 y.n+1/2 j^g^g^jj ^^.^^ ^^lat the physical flux F{W) function for the 
acoustic equation (2.4.5) is equal to 

(2,4.19) F(W) = (^Pf;) with W= 

Proposition 2.5. Upwind scheme for computational acoustics. 

The extension of the upwind finite volume scheme (2.3.10), (2.3.12) and (2.3.14) 
is determined by the following relations : 

r2 4 20^ .n+l/2 ^ / ^ (wy_i/2 + - f{P'j + l/2-P]-l/2) 

\ {Pj-l/2+Pj+l/2) - -^j-1/2) 

for the internal fluxes, i.e. for indexes j that satisfy 1 < j < J—1. The two 
boundary fluxes follow the following relations : 

.n+1/2 ^ f + Co (n((n+ i)At) -P1/2) 

V ^ in((n+i)At) 

^2 4 22"! fri+1/2 _ I (wj_i/2 + ^j-1/2) ~ ^ (P^J -1/2 ~ P"} -1/2) 

\ \Pj-l/2 '^Pj-1/2) ~ 2 l^J-1/2 ~ ^J-1/2/ 

• The internal fluxes are determined with the scheme (2.3.10) applied with 
the diagonal form of relation (2.4.9). We have 

(2.4.23) = = P"-i/2 + Poco 

(2.4.24) v-tf ' = V^,,+i/2 = P"+i/2 - Poco 
then the relation (2.4.20) is established. 

The left boundary flux uses the extension of relation (2.3.12). We first determine 
the characteristic variables on the left boundary according to relation (2.4.13) 

(2.4.25) vl-^i'' = 2n((n + i)At) - <^:«/^ 

and use a first order extrapolation of the outgoing characteristic variable : 

(2.4.26) ipV:i^' = ipl^y^ = K/2 - Poco </2. 

Then we solve the system (2.4.25) (2.4.26) and find finally the relation (2.4.21). 
The process is analogous for the right boundary. The input datum is imposed 
according to the relation (2.4.14) : 

(2.4.27) ^l~iy^ = = Po{L) - poCoUo{L) « P^j-1/2 - PoCq u^j_i/2 
and the output characteristic variable is extrapolated from the interior of the 
domain : 

(2.4.28) (^;+y' = (^^,j_i/2 = p}-i/2 + Poco n}_i/2. 

The relation (2.4.22) follows after two steps of elementary algebra. □ 
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• We remark that both relations (2.4.20) and (2.4.22) are identical, except 
that the boundary state Wo{L) ^ ^las replaced the right state W'^_^^i2 ■ 

Moreover the flux boundary condition (2.4.21) that involves the pressure is a 
natural discretization of the exact characteristic solution (2.4.17) (2.4.18) at 
x = {). 



2.5 Characteristic variables. 

• We suppose now to flx the ideas that the unknown vector W(«, •) 
(2.5.1) [0, L] X [0, +oo[ 3 (x-, t) I — y W{x, t) e 

has three real components wi, W2 and ws. We suppose also that the function 
W{», •) is solution of a conservation law of the type 
dW d 

where the flux F{W) is a linear function of vector W : 

(2.5.3) F{W) = A.W 

and A is a 3 by 3 diagonalizable real matrix. 

• We flrst detail the fact that matrix ^ is a diagonalizable matrix. There 
exists three non null real vectors ri , r2 , rs and three real scalars Ai , A2 , 
A3 in such a way that 

(2.5.4) A.rj = X,r, , j = 1, 2, 3. 

^Prom a matricial viewpoint, we denote by Rkj the component of the 
eigenvector rj , i.e. 

(Rij\ 




(2.5.5) r 

and we introduce the 3 by 3 matrix R composed by the scalars Rkj ■ The vector 
Tj is the column of matrix R. The relation (2.5.4) can also be written as 

(2.5.6) A.R = R.A, 

and A is the diagonal matrix whose diagonal terms are equal to the eigenvalues 
A,- : 

/Ai 

(2.5.7) A = A2 

V A3 

• We consider now two distinct bases for linear space IR"^ : on one hand 

the canonical basis (e^) ^ ^ deflned by 

/1\ ' ' /O 

(2.5.8) ei = , 62 = 1 1 , 63 
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where the vector W admits the natural decomposition introduced above : 

(2.5.9) W = YiZi yokek , 

and on the other hand the basis of IR^ composed by the eigenvectors {rj)j=i^ 2, 3- 
In the latter, the vector W can be decomposed with a formula of the type 

(2.5.10) W = E^l^jT, 

and the scalar (pj define the characteristic variables associated with the 
system (2.5.2) (2.5.3). The link between the relations (2.5.9) and (2.5.10) is 
classical : we consider the components R^j of vector rj inside the canonical 
basis and we get from the relation (2.5.5) : 

(2.5.11) Wk = YljZl^jRkj- 

Then the relation (2.5.11) can be re- written under a matricial form : 

(2.5.12) W = R.if. 

• The relation (2.5.12) proposes to change the unknown function, i.e. 

to replace the research of W{x^t) e IR'^ by the equivalent research of the 
characteristic vector (/?(a:, t) G IR^ and defined by : 

(2.5.13) ^ = R-^.W. 

Proposition 2.6. Characteristic variables satisfy advection equations. 

The vector [0, L] x [0, +00 [ 3 (x, t) 1 — y ip{x, t) G IR"^ of characteristic variables 
satisfy the matrix equation 

(2.5.14) t.A.g^O 

that takes also the equivalent scalar form : 

(2.5.15) ^ + A,^ = 0, i = l,2,3. 



• We have from (2.5.2), (2.5.3), (2.5.6) and (2.5.12) : 

+ A—- = R.^ + A.R.^ = R.l^ + R-^.A.R.^ 

at ox ot ox \ ot ox J 

and since the matrix R is invertible, we deduce from the previous calculus the 
relation (2.5.14). The relation (2.5.15) is an immediate consequence of (2.5.14) 
and (2.5.7). □ 
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/ 



X 3 



/ ^ 



x = 



X 3 



/ X 



x = L 



Figure 2.7. Linear hyperbolic system with three equations 
and eigenvalues satisfying Ai < < A2 < A3. 

• To fix the ideas, we suppose that the eigenvalues Aj of matrix A 
are distinct, enumerated with an increasing order and with distinct signs as 
illustrated on Figure 2.7 : 

(2.5.16) Ai < < A2 < A3. 

The propagation of the first variable goes from right to left (because Ai < ) 
with celerity | Ai |, the second characteristic variable (/?2 from left to right with 
celerity A2 and the same property holds for variable (/?3 with eigenvalue A3. 

• A set of well posed boundary conditions is a consequence of the diagonal 
form (2.5.15) of the equations and of the particular choice (2.5.16) for the signs. 
The directions associated with eigenvalues A2 and A3 are ingoing at a; = and 
we have to give some boundary condition for (/?2 and (/?3 at this point : 

(2.5.17) ¥^2(2^ = 0, t) = /5o(t) 

(2.5.18) ¥^3(a: = 0, i) = 7o(t). 

The direction associated with the eigenvalue Ai is ingoing at the abscisssa x = L, 
and this condition imposes to have some datum concerning (pi at this particular 
point : 

(2.5.19) ipi{x = L, t) = aL{t) . 

The previous boundary conditions (2.5.17) to (2.5.19) define a well posed prob- 
lem. Nevertheless, the introduction of physically relevant boundary conditions 
(as a pressure condition as seen in the previous section) requires a more general 
formulation of the boundary condition. In the linear case, the stability study 
developed by Kreiss [Kr70] shows that the ingoing characteristic can be an affine 
function of the outgoing characteristic through a reflection operator at the 
boundary. We can explicit the former with the above example. 
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Figure 2.8. Reflection operator at x = 0. 




Figure 2.9. Reflection operator at x = L. 

• At a; = , the first characteristic is outgoing and the two last ones are 
going inside the domain of study. Then we can replace the conditions (2.5.17) 
and (2.5.18) by the following ones : 

(2.5.20) if2{x = Q>,t) = Po{t) + p{t)ipi{x = 0,t) 

(2.5.21) ip3{x = 0,t) = 7o(t) + q{t)ipi{x = 0,t) , 

where t i — > p{t) and t i — > q{t) are given fixed real functions of time. The 
conditions (2.5.20) and (2.5.21) are illustrated on Figure 2.8. We can also write 
them 

(2.5.22) (fi'^^ix^t) = g{t) + S (t) • (p''''* {x , t) , a; point on the boundary, 

with v.-(::)-(*)-(5g)'^<*^^(S))-°"'^- 

when X = . 

• When x = L, the relation (2.5.19) is replaced by a more general one 

(2.5.23) ipi{x = L, t) = aL{t) + e{t) ^2{x = L, t) + (j(t) ^z{x = L, t) 
illustrated on Figure 2.9 and including an affine component of the outgoing 
characteristic variables. The boundary condition (2.5.23) takes again a form of 
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the type (2.5.22) with this time the following relations : 99*" = (pi , g{t) 
aL[t) , S{t) = (e{t) a{t)) , if"''^ = f when x = L. 



2.6 A family of model systems with three equations 

• We still study a 3 by 3 linear hyperbolic system of the type (2.5.2) 

(2.5.3) with the condition (2.5.16) to fix a particular example. We suggest in 

this section to explicit a way for evaluation of the numerical flux f^^^^'^ that 
is the key point for the discrete evolution in time of the mean values Wj+1/2 : 

(2.6.1) ^(1^1% -wj^y.) + ^^{fj:^" - ir") = 0- 

The internal fluxes {f^~^^^'^) are evaluated with the help of a two- 

point numerical flux function $(•, •) : 

(2.6.2) /;+^/^ = $(w7_i/2, w-;+i/2) 

undary fluxes J 

sub-section 



and the boundary fluxes Jq'^^^'^ and are detailed in a forthcoming 



t 



w 

left 



right 



T 

Figure 2.10. Discontinuity at the interface between two cells. 



• We change the notations and wish to determine the numerical flux 

^{Wi, Wr) for Wi = Wieft and Wr = Wright given respectively at the left 
and at the right of the interface (see Figure 2.10). When we consider the ad- 
vection equation (and in that case the variables Wi and Wr are real num- 
bers) the relation (2.1.19) gives the result : ^{Wi, Wr) = aWi when a > 
and ^{Wi, Wr) = aWr when a < 0. We have to generalize this study when 
the field W{», •) is three-dimensional. We first decompose the vector $(W/, 
Wr) with the basis rj of eigenvectors and introduce its (scalar) components 
^j{Wi, Wr) : 

(2.6.3) ^Wl, Wr) = E^=? MWl, Wr) Tj 

i.e. 

(2.6.4) ^k{Wi,Wr) = E^lRkjiJj{Wi,Wr). 

For J = 1, we have Ai < then the numerical scheme has to be upwinded in 
the right direction : 
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(2.6.5) ^iW, Wr) = \l^l,r 

whereas for j = 2 or ^' = 3, we have A2 > and A3 > and the scheme must 
be upwinded to the left. It comes 

(2.6.6) ij2{WuWr) = \2^2,l. MWuWr) = Xs^S,!- 

In consequence of the relations (2.6.3) to (2.6.6), the numerical flux function 
$(•, •) can be written globally : 

(2.6.7) ^{Wi,Wr) = Xi(fi,rri + X2(p2,ir2 + X3(f3,ir3 , 

or in an equivalent way with introducing the Cartesian components : 

(2.6.8) ^k{Wi,Wr) = Xiipi^rRkl + X2ip2,lRk2 + X3ip3,lRk3, A;=l, 2, 3. 

• We can also re-write the relation (2.6.8) for the particular interface Xj : 

(2.6.9) Wi = W^ieft = W7_i/2 , Wr = M^right = . 

We first decompose the vector W on the eigenvectors of matrix A as in (2.5.11) : 

(2.6.10) (^;+i/2)fc = E-3^.^,+i/2^fci, A;=l,2,3, j = 1, • • • , J-1 , 
then we introduce the component number of the flux /"^^^^ , i.e. {f^'^^^^)k 
= ^fc(W'j^_j/2' ^j+1/2) interface xj : 

(2.6.11) (/r^"^^^)fc = ^1 ^1, j + 1/2 ^fcl + A2 (^2, j-1/2 -^fc2 + A3 (^3' j_l/2 ^A;3 • 

• We detail in this sub-section the determination of the numerical flux 
jn-i-1/2 boundary x = 0. We first recall that the continuous boundary 
conditions at this point take the form given in (2.5.20) (2.5.21). The idea is to try 
to apply the upwind scheme (2.6.11) at the particular vertex ^' = : fg^^^"^ = 
Ai (/?^ ^1 + + -^2 ^2 -1/2 ^2 + A3 (yC'g _^^2 r3 and then to replace the charac- 
teristic values _i/2 ^3 -1/2 (that are not defined on the mesh) by their 
values evaluated after a rough discretization of relations (2.5.20) and (2.5.21) : 

_ nn+l/2 , n+1/2 n n _ ri+1/2 n+1/2 n 

^2,-1/2 — PO ^ '^1,1/2' '^3,-1/2 —10 ^ y '^1,1/2- 

We obtain in consequence the following expression for the boundary flux at 
x = : 

,n+l/2 _ / Al Vl 1/2 n + A2 ( /3r^/' + Vl 1/2 ) '■2 + 

• ' ^° "I +A3(7o""'^^ + 9"+'/V?,i/2)r3 
or in an equivalent way : 

(2 6 13) r+^/' = / V2 (^1^1 + ^2 P"-^'^' r2 + A3 r3 ) + 

^ ■ ■ ^ ' 1 +A2/3?+^/^2 + A3 7?+'/'r3. 

• The determination of the boundary flux can be conducted 
in the same way. Starting from the expression of the upwind scheme (2.6.11) 
when j = J, i.e. formally jy+V2 ^ Xiip^^j^^j^ri + A2 (^2, j-1/2 ^2 + 



26 



An introduction to finite volumes for gas dynamics 



^3 <^3 J -1/2 ^3 ' replace the first characteristic variable that appears ex- 
ternal of the domain by its value given by the boundary condition (2.5.23) : 

Kj+1/2 = + ^"+'/V^,j_i/2+ +^"+'/V^,j_i/2- We deduce : 

(2 6 14) = / ^1 ( + ^2, j-i/2 + j-i/2 ) ri 

\ +M^2, J-1/2 ^2 + A3 ifl rs 

or in an equivalent manner : 

(2 6 15) r^''^ = + ^2, j-1/2 ( Ai n + A2 r2 ) + 

1 +^3,j-i/2(Aia"+V'ri + Asrs) . 

2.7 First order upwind- centered finite volumes 

• Wc consider now a general system of conservation laws 

dW d 

with an unknown vector •) that belongs to linear space IR"^ : 

(2.7.2) [0, L] X [0, +oo[ 3 {x, t) ^ W{x, t) G IR^ 
and a linear flux function F{») 

(2.7.3) F{W) = A.W 

associated with a diagonalizable matrix A with eigenvalues \j and eigenvectors 



^3 



(2.7.4) A.Tj = XjTj , j = 1, 2, m. 

Introducing the mxm matrix R as in relation (2.5.5) and the diagonal matrix 
A of eigenvalues as in relation (2.5.7), we have : 

(2.7.5) A.R = R.A. 

• We propose here to determine a first order upwind flux ^(Wi,Wr) be- 

tween the two states VFieft = Wi and Wright = Wr that generalizes the relation 
(2.6.7) when we have not done any hypothesis of the type (2.5.16) concerning 
the sign of the eigenvalues Xj. We decompose any state W on the basis of space 
IR"^ characterized by the eigenvectors rj : 

(2.7.6) W = Y.T=i ^jTj , Wi = ^3,irj , Wr = E^ZT 



j=l "TJ • 3 1 I — Z^j=l "^3,1' 3 1 ""r ~ ^j=l "t^3,'r'^3'> 

and due to the structure introduced at Proposition 2.6, we obtain an advection 
equation for the j° characteristic variable : 

Therefore it is natural to introduce the components ipjiWi, Wj) of the numer- 
ical flux on the basis of the eigenvectors : 

(2.7.8) ^{Wi, Wr) = EjS V^^W, Wr)rj 
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and the first order upwind finite volume scheme is defined by the way we eval- 
uate the coefficient ipjCWi, Wr) with the upwind scheme associated with the 
advection equation (2.7.7) : 

(2.7.9) ^,(WuW^ = l^'^;^ ^^11 

• For any real number fi , we introduce the positive part and the 
negative part fj,~ by the relations 

^' = {'o if^fo' ^- = {1 if^i:. 

We remark that we have 

(2.7.11) = , V/xGlR 

(2.7.12) =//+-//- , V//GIR. 

We introduce also the absolute value | A | of the diagonal matrix A by the 
condition : 

(2.7.13) I A| = |diag(Ai, • • • , Am) I = diag( | Ai |, • • • , | | ) 

and due to the relation (2.7.5), the absolute value \ A \ of the matrix A is 
defined by : 

(2.7.14) 1^1 = R. |A| .R-^ . 

Proposition 2.7. Three expressions of the upwind first order scheme. 

Let ^{Wi, Wr) the upwind flux defined by the relations (2.7.8) and (2.7.9). 
Then we have the three equivalent expressions : 

(2.7.15) ^Wi, Wr) = F{Wi) + ^7 (^j,r - rj 

(2.7.16) ^Wi, Wr) = F{Wr) - TizT A+ (^j,r - rj 

(2.7.17) ^Wi, Wr) = \ (FiyVi) + F{Wr)) - \ \A\ .{Wr - Wi) . 

• We write the relation (2.7.9) under the form : 

(2.7.18) i),{Wu Wr) = A+ ipj,i + AT ipj^r 

and we have : 

^W, Wr) = (K ^j^l + V '^j^r)rj 

= Ej=r [i^j - \") + V V^j,r)rj due to (2.7.11) 

and the relation (2.7.15) is established. In an analogous way, we have : 
^{Wu Wr) = E'jZT (A+ + AT ^,,r)r, 

= Ej=r + (A, - A+) ip,,r) r, due to (2.7.11) 

^{Wu Wr) = T.i=T A, ^,,r T, - Y.]Z7 A+ {^j,r ' ^j,l)r, 

and the relation (2.7.16) holds. We remark that 
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\A\ •{Wr-Wi) = R* |A| .R-^ .R.{ifr-^i) due to (2.7.14) and (2.5.12) 

= R* I A| •{Lfr - (fl) 

= Ejir ^fcj l^jl i^j,r-<fj,i)ek then 

(2.7.19) \A\.{Wr-Wi) = I A,- 1 {^j,r-^j,i)rj. 

We add the previous results (2.7.15) with (2.5.16), and we divide by two. We 
obtain : 

^Wl, Wr) = I (F{Wl) + F{Wr)) - \ Y.T=1 (^t - V) (^j^r - Tj 

= \ {F{Wi) + F{Wr)) - \ Ej=r I A,- 1 fe,. - Tj due to (2.7.12) 
= \ + F{Wr)) - \ \A\ .{Wr - Wl) 

due to the relation (2.7.19). Then the relation (2.7.17) is established and the 
proposition 2.7 is proven. □ 



3) Gas dynamics with the Roe method. 

3.1 Nonlinear acoustics in one space dimension. 

• We propose here to describe quickly a physical problem that comes 
from the theoretical modelling of trombone, detailed for instance in the work 
of Hirschbcrg et al [HGMW96] or in our study [MD99] with R. Msaham. In a 
first approximation, the duct of a trombone is a long cylinder with a constant 
section and the acoustic waves propagate only in the longitudinal direction. We 
can use a one- dimensional description of the geometry (see Figure 3.1) and in 
what follows, the trombone is modelled by a real space variable x that ranges 
from x = Q at the input to x = L at the output. 

• At the input a; = 0, a given non-stationary pressure wave ti — >■ n(t) is 
emitted ; this wave is a perturbation of the ambiant pressure po of the air : 
(3.1.1) |n(t)-2?o| « , t>0. 

At the output X — the waves go outside without any reflection due to the 
presence of a pavilion and the boundary condition is a "free output" and a non- 
reflecting boundary condition has to be used. At the initial time t = 0, we 
can consider that the air satisfies the usual conditions of pressure p{x, 0) = po , 
temperature T{x^ 0) = To and density p{x, 0) = po- We study in this section a 
finite volume method able to treat nonlinearities in the acoustic modelling and 
based on the characteristic decompositions developed in the previous section. 
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Input nonreflecting 

pressure free 

n(t) output 
x=0 x=L 

Figure 3.1. Long unidimensional pipe for the modelling of a trombone. 



3.2 Linearization of the gas dynamics equations. 

• We study a perfect gas subjected to a motion with variable velocity in 
space and time. We have noticed that the primitive unknowns of this problem are 
the scalar fields that characterize the thermodynamics of the gas, i.e. density p, 
internal energy e, temperature T, and pressure p. In what follows, we suppose 
that the gas is a polytropic perfect gas ; it has constant specific heats at constant 
volume Cv and at constant pressure Cp. These two quantities do not depend 
on any thermodynamic variable like temperature or pressure ; we denote by 7 
their ratio : 

C 

(3.2.1) 7 = (= constant) . 

We suppose that the gas satisfies the law of perfect gas that can be written with 
the following form : 

(3.2.2) p = (7-l)pe. 

As usual, internal energy and temperature are linked together by the Joule- 
Thomson relation : 

(3.2.3) e = CyT. 

• In the formalism proposed by Euler during the IS*'^ century, the motion 
is described with the help of an unknown vector field u which is a function of 
space X and time t : 

(3.2.4) u = u{x, t) . 

In the following, we will suppose that space x has only one dimemsion [x £ 
IR). We have four unknown functions (density, velocity, pressure and internal 
energy) linked together by the state law (3.2.2). In consequence, we need three 
complementary equations in order to define a unique solution of the problem. 
The general laws of Physics assume that mass, momentum and total energy are 
conserved quantities, at least in the context of classical physics associated to 
the paradigm of invariance for the Galileo group of space-time transformations 
(see e.g. Landau and Lifchitz [LL54]). When we write the conservation of mass, 
momentum and energy inside an infinitesimal volume dx advected with celerity 
w(a;, t), which is exactly the mean velocity of particules that compose the gas, 
it is classical [LL54] to write the fundamental conservation laws of Physics with 
the help of divergence operators : 
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(3.2.6) + ±{pu^+p) = 

• We introduce the specific total energy E by unity of volume 

(3.2.8) E = ^u"^ + e , 

the sound celerity c following the classical expression : 

(3.2.9) c = 

and total enthalpy H defined according to 

(3.2.10) H ^ E + ^^ = ^u^ + ^^c\ 

The vector W is therefore composed by the "conservative variables" or more 
precisely by the "conserved variables" : 

(3.2.11) W = (p,pu,pE)^ = {p,q,e)\ 

The conservation laws (3.2.5)-(3.2.7) take the following general form of a system 

of conservation laws : 

dW d 

where the flux vector W 1 — > F(W) satisfles the following algebraic expression : 

(3.2.13) F{W) = {pu, pu^ + p, puH)^ 

that can be explicited as a true function of state vector W, on one hand with 
the pressure law P{W) computed with (3.2.2), (3.2.8) and (3.2.11) : 

(3.2.14) P{W) = {^-l)(^e-Q 

and on the other hand with an explicit use of the conserved variables p, q and 
e. We obtain : 

2 

(3.2.15) F{W) = (q,^+PiW),^+P{W)^). 

\ P p py 

Proposition 3.1. Jacobian matrix of gas dynamics. 

• The Jacobian matrix dF{W) of the flux function W 1 — > F{W) for the 
Euler equations of the gas dynamics admits the following expression : 

/ 10 

(3.2.16) dF(^) = (7-l)i/-w2-c2 (3 - 7) w 7-I 

\ {-f-2)uH -uc^ H-{j-l)u'^ -fu 

• The matrix dF(W) is diagonalizable ; the eigenvalues Xj{W) satisfy the 
relations 
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(3.2.17) Ai(T^) = M-c < X2iW) = u < Xs{W) = u + c. 

and the associated eigenvectors rj{W) are proportional to the following ones 

f ' \ ( '\ ( ' 

(3.2.18) ri{W)= u-c , r2{W) = \ u \ , r^{W) = u + c 

\H-uc) \h + uc 

• We first differentiate the pressure law W i — )■ P{W) given in (3.2.14) : 

, dP 7-1 9 , , n dP , , dP , , 

^^•^•^^^ a;^ ^ 2 ^ ^ ' ^ = ^ = 

and the second row of the matrix (3.2.16) is a direct consequence of the relations 

• The calculus of the third row of matrix in (3.2.16) demands first evalu- 
ation of the gradient of puE = ue relatively to the state W. We get 

(3...0) |(^)^-«., |(^)^., |(^)^.. 

We have also -^{Pu) = + p-^{^) then we deduce from (3.2.19) 

and the following expressions for the gradient of velocity ^(p) = and 

d /'Pq\ 7 — 1 q up 



dq\pj p 



(3.2.21) { ^riS ^ ^ P 

We add the relations (3.2.20) and (3.2.21) ; then the third row of matrix (3.2.16) 
admits the following expression : ( ^ — uH ^ iJ - (7— l)w^, 7^) and 
this result is exactly the third row of the right hand side of (3.2.16) when we 
take into account the relation (3.2.10) between and c^. The relations 

(3.2.17) and (3.2.18) are elementary to satisfy ; they express simply the three 
relations : 

(3.2.22) dF{W).rj{W) = \j{W)rj{W) , j = 1, 2, 3 

and Proposition 3.1 is established. □ 

• We keep into memory the following expression of the Jacobian matrix 

dF{W) : 

/ 10 

(3.2.23) dF{W) = (3-7)w 7-I 

\^u^ -uH i?- (7-1)^2 

that needs only the datum of velocity u and total enthalpy H of the state W. 
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3.3 Roe matrix. 

• We consider two states W\eit = Wi and T^right = Wr relatively to the 

gas dynamics, i.e. they both belong to space JR'^ and have an expression of 
the form (3.2.11). By definition, a Roe matrix A{Wi, Wr) between these two 
states is a 3 by 3 matrix that satisfy the three following properties : 

(3.3.1) A{Wi, Wr) is a diagonalizable matrix on the field IR of real numbers 

(3.3.2) A{W, W) = dF{W) 

(3.3.3) F{Wr) - F{Wi) = A{Wi, Wr) . {Wr - Wi) . 

In his original article, P. Roe [RoeSl] has proposed a very simple algebraic way 
to construct a Roe matrix for the dynamics of polytropic gas. We propose it in 
the following Proposition. 

Proposition 3.2. Algebraic construction of a Roe matrix [RoeSl] . 

Let Wi and Wr be two states for gas dynamics, defined by their densities pi 
and Pr, their velocities ui and Ur and their total enthalpies Hi and Hr- We 
introduce an intermediate state VF*(VF;, Wr) by its density p*, its velocity 
u* and its total enthalpy H* according to the following relations : 

(3.3.4) p* = ^ypl pr 

(3.3.5) u = 



(3.3.6) i?* = ^ ^ - 

Then the matrix A(W/, Wr) defined as the Jacobian matrix of the flux for the 
intermediate state W*(VFz, Wr), i.e. 

(3.3.7) A{Wu Wr) = dF(W*{Wi, Wr)) 
is a Roe matrix. 

• Due to the expression (3.2.23) of the Jacobian matrix of gas dynamics, 

we remark that the formula (3.3.4) giving the density p* is not necessary for the 
determination of the matrix dF{W*{Wi, Wr)) and an entire family of states 
VF*(VF/, Wr) define a Roe matrix according to the relations (3.3.5), (3.3.6) and 

(3.3.7) . Nevertheless, we keep this definition of density p* by convenience and 
simplicity for future algebraic expressions. The proof of Proposition 3.2 needs 
some algebraic developments. We begin by the following technical lemma. 

Proposition 3.3. 

Under the hypotheses of Proposition 3.2, we have the following relations : 

(3.3.8) {U*)'^{pr- pi) - 2U* {prUr- plUi) + {prUl- piuf) = 

(3 3 9) I ^* ~ ^ ^* ~ ^ ^* ~ ^ 

1 = prUrHr - plUiHi . 
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• We first evaluate the left hand side of relation (3.3.8) : 

{U*f [pr- pi) - 2 U* {pr Ur- plUl) + {pr - pi uf) = 

= U*{^/p^-^/p'i) {^fp^Ur + ^/plUi) - 2u*{prUr-plUi) + {prul-piuf) 



U 



iVpi {VPI + VP^) Ul - yfp^ {yfpl + ^fp^) Ur) + (pr 



Pi Ul) 



= {^lUl + y/p^Ur) {^/plUl - ^rUr) + {pr - PlUi) 

= and the relation (3.3.8) is established. 

We work on the left hand side of (3.3.9) as follows : 

-U* H* {pr - pl) + H* {pr Ur - pi Ui) + U* {pr Hr - pi Hi) = 

-- - U* {^f p^ - ^i) {y/plHi -\-y/p^Hr) + H* {pr Ur - pl Ui) + U* {pr Hr - pl Hi) 
-- VPI Pr U* {Hr - Hi) + H* {pr Ur - pl Ui) 

\fpi^J~P'r{^fPlUl + yfp^Ur) {-Hi +Hr) + {- pl Uj + Pr Ur) {y/pl Hi + ^Hr) 

\fpl+ ^ 

1 

[ -pl {VPI + VP^) ^l Hi + pr iVPl + y/p^) Ur Hr ] 



y/pi+VP^ 
= Pr Ur Hr - pl Ul Hi 

and the proposition 3.3 is established. 



□ 



• The proof of Proposition 3.2 consists in satisfying the three hypothe- 
ses that define a Roe matrix. First, due to the fact that the relation (3.3.7) 
defines the matrix A{Wi^ Wr) as a Jacobian of some state, this matrix is diag- 
onalizable with real elements due to the result of Proposition 3.1 and the first 
property (3.3.1) is satisfied. The second property (3.3.2) is a simple consequence 
of the fact that if Wi = Wr = VF, then we have from the relations (3.3.4) to 
(3.3.6) : W*{Wi, Wr) = W and the property results from (3.3.7). 

• The third property (3.3.3) needs more work. We remark that the first 
row of this matricial relation is clear. For the second row, we have : 

Second row of matrix A{Wi, Wr) • {Wr -Wi) = 



7-3 



{U*)'^ {Pr - pl) + (3-7) M* {prUr-plUi) + {'y -1) {pr Er - pl Ei) 



2 

[(W*)^ {pr-pl) - 2 U* {pr Ur - pl Ul)] H — {pr U^ 



= {pr ul - pi uf) + {pr - Pl) 

= second row of the flux difference F{Wr) — F{Wi) 

• We have also, in consequence of (3.2.23), 
Third row of matrix A{Wi, Wr) • (Wr -Wi) = 
-7-1 



PlUi) + {Pr-Pl) 

due to (3.3.8) 



u (- 



7-1 



(W*)2 - H*) {pr -pi) + {H*- (7-1) (n*)2) {pr Ur - Pl Ul) + 
+ 7W* {prEr-piEi) 



u 



U*)'^ {pr-pl) - 2 U* {pr Ur - pl Ui) \ + 
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+ [-H* U* {pr - Pi) + H* {pr Ur - pi Ui)] + 7 {pr Er - pi Ei) 

7— 1 

= — ^ U* [{U*Y [pr - pi) - 1 U* [pr Ur - pi Ul)] - U* {pr Hr - pi Hi) + 

+ {pr Ur Hr-pi Ui Hi) + 7 u* {pr Er-pi Ei) due to (3.3.9) 

= U* [ (U*)^ {pr- pi) - 2 U* {pr Ur- plUl) + pr U^ - pi uf ] + 

+ W* {-JPrCr +^ piei +^ prCr - J Pi Ci) + prUrHr - plUiHi 

= prUrHr—piUiHi duc to (3.3.8) 

= third row of the flux difference F{Wr) — F{Wi) 

in the view of relation (3.2.13). The proposition 3.2 is established. □ 
3.4 Roe flux. 

• The principal interest of the Roe matrix is to be able to use all what has 

been developed for linear hyperbolic systems in Section 2. In particular, the 

following linear hyperbolic system defined with a given Roe matrix A(W/, Wr) 

dW dW 
(3.4.1) — + A(W„W^).-^, 

can be treated with the upwind scheme defined at proposition 2.7. We obtain 
by doing this the following 

Proposition 3.4. Three formulae for a flux. 

• Let Wi and Wr be two fluid states and W* the intermediate state 
defined by the relations (3.3.4) to (3.3.6). The sound celerity c* of state W* 
is defined with the help of relation (3.2.10), i.e. 



(3.4.2) c* = Y^(7-l)(i?*-^^ 



and the eigenvalues A* of the Roe matrix A{Wu Wr) = dF{W'{Wi, Wr)) 
are given by a relation analogous to (3.2.17). 

(3.4.3) Xl=u* -c* < X^ = u* < = + c* . 

The associated eigenvectors = rj{W*) are proportional to the following 
ones : 

/ 1 \ /I 

(3.4.4) rl = 



u*-c* , 



\H* -u* c* 



u* , r* 



VIM' 



3 




We introduce the decomposition of vector Wr — Wi in the basis r*- 



(3.4.5) Wr-Wi = E^=? «i . 

The three following relations define a unique numerical flux $(W/, Wr) named 
the Roe flux between the two states Wi and Wr : 

(3.4.6) ^Wi, Wr) = F{Wi) + Y^^Zl (A*)- aj r* 

(3.4.7) ^Wi, Wr) = F{Wr) - Ej=? (A*)+ cij r* 
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(3.4.8) HWi,Wr) = ^ (F{Wi) + F{Wr)) - ^ \A{Wi,Wr) \ .(Wr-Wi). 
• The first non-obvious point is to verify that the relation (3.4.2) defines 



a real number c*. We have 



*\2 



2 




1 








1 






/p-r? 


1 





^iHi -\- y/p^Hr I ( ^lUi + ^/p^Ur \^ 

2 



{^/pi+Vp~r) (v^(2^?+:^^n + Vp^ {2^1+7— ["^D) 

- ^ (P« + 2 P* Ui Ur + pr 



{\fPl + xfP^)'' 



\[p* U^l -2p* UiUr + p* ul) + 

1 * / x2 , + 2 , P* + pr 2 
- /? {Ur- UiY + ^ Cf H — c: 



P/+P* 2 , P*+P 



r 2 



7-1 



7-1 



7-1 
> 



and 
(3.4.9) 



c* = 



^^p*{Ur-UiY + {pi + p*) + {p* + Pr) 



\fpi + \/P^ 

• We make the difference between the right hand sides of (3.4.6) and (3.4.7). 

We get : 

F{yVr)-F{yVi) - E^C?((A*)+ + (A*)-) r* = 

= AiyVu Wr) . iyVr - Wi) - Y^^r^x a* aj r* due to (3.3.3) and (2.7.11) 

= Aiyvu Wr) . (Ej3 «i - Y.T=-y ^* ^1 (3.4.5) 

= because AiyVi^ WV)»r^ = A* r * for each integer j. 

The proof of relation (3.4.8) is obtained by taking the half sum of (3.4.6) and 
(3.4.7). It is analogous to the one done for Proposition 2.7. The proof of Propo- 
sition 3.4 is completed. □ 

• We make explicit the parameters aj introduced in relation (3.4.5) 
in order to be complete for the implementation of the above formulae on a 
computer. 

Proposition 3.5. New acoustic impedance. 

With the notations introduced at Proposition 3.4, and denoting by and Pr 
the respective pressures of states Wi and Wr-, we have the following relations for 
the scalar components aj of the state difference Wr — Wi in relation (3.4.5) : 



(3.4.10) 
(3.4.11) 



Q!2 



2(c*)2 

1 



[ {Pr - P* C* Ur) - {pi - p* C* Ul) ] 



iPr - {cn' Pr) - {Pl - {cn' Pl)] 
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(3.4.12) as = [ {Pr + P* c* Ur) - {pi + p* c* Ul) ] 

with acoustic impedance p* c* that is nomore the one po Cq of a reference 
state as in traditional acoustics but an impedance associated with the Roe in- 
termediate state W*{Wi,Wr) of relations (3.3.4) to (3.3.6). 

• We have just to explicit the three components of the relation (3.4.5). It 
comes : 

Pr- Pl \ / ttl + ^2 + «3 

PrUr-plUl = CKi (W* -C*) + a2W*^+ a3(w* + C*) 

PrEr-piEi) \ai{H* -u*c*) + + as{H*+u*c*) 

(3.4.13) ai+a2 + as = Pr - Pi , 

and we deduce after multiplying the equation (3.4.13) by —u* and adding to 
the second equation of the above matrix equality : 

C* (^3 - ai) = prUr - pl Ui - U* [pr - pl) 

= prUr - pl Ul - {^fp^ - ^/pl) {^lUi + ^/p^Ur) , 

then 

(3.4.14) C* (— tti + tts) = p* {Ur — Ul) . 

• We deduce from the third equation of relation (3.4.5) : 
(c*)2 1 

{ai + as) = prEr - piEi - - (m*)^ {pr - Pl) - u* c* {as - ai) 

'f — L Z 

= TiPr-Pl) + -(PrUl-piuf) + 7:K)^(Pr-Pz) - U* {pr Ur - pl Ul) 

■7—1 z z 

= {pr — Pl) due to (3.3.8). Then we have : 

7-1 

(3.4.15) (c*)2(ai+a3) = Pr-Pi. 

• The solution of the 2 by 2 linear system with unknowns ai and as 
defined by the relations (3.4.14) and (3.4.15) directly gives the relations (3.4.10) 
and (3.4.12). The expression (3.4.11) of variable a2 is a direct consequence of 
the relations (3.4.10), (3.4.12) and (3.4.13) and Proposition 3.5 is proven. □ 

Proposition 3.6. An algorithm for the Roe flux. 

Let Wi and Wr be two compressible fluid states. The computation of the 
Roe flux ^{Wu Wr) of relations (3.4.6)-(3.4.8) between these two states is 
summarized by the following points : 

• Evaluation of density p*, velocity u* and total enthalpy H* of the 
intermediate state W* with the relations (3.3.4) to (3.3.6), 

• Determination of the sound celerity c* of the intermediate state W* from 
the previous data with the relation (3.4.2), 

• Eigenvectors r*j of the Roe matrix with the relations (3.4.4), 
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• Computation of the characteristic variables aj in (3.4.5) for the difference 
Wr - Wi with the relations (3.4.10) to (3.4.12), 

• Final computation of the Roe flux $(W/, Wr) with the minimum of work : 

(3.4.16) ^{Wi, Wr) = F{Wi) if w* - c* > 

(3.4.17) $(W;, Wr) = F{Wi) + (w* - c*) tti if u* - c* < < u* 

(3.4.18) $(Vz, Wr) = F{Wr) - (u*+c*) as if < < n*+c* 

(3.4.19) ^{Wi,Wr) = F{Wr) ifi/* + c*<0. 

• The proof of the relations (3.4. 16)- (3.4. 19) is obtained by starting from 
the expression of the Roe flux given in (3.4.6). We know that Ai = w* — c*, 
A2 = u*, Xs=u* + c*. If u* -c* > 0, then u* >0 and w* + c* > 0, so the 
relation (3.4.6) reduces to (3.4.16) because, due to (2.7.10), = if // is a 
positive real number. If u*—c*<0 < u* <u*-\-c*, the term containing (A^)~ 
is the only one that contributes in relation (3.4.6) and the relation (3.4.17) is a 
direct consequence of this remark. If u* — c* < u* < < + c*, the term 
that contains (A3)+ is the only nonzero element among the three inside the 
relation (3.4.7) and we deduce the relation (3.4.18) from this property. When 
u* -c* < u* < u* +c* < 0, the term F{Wr) is the only to subsist inside the 
relation (3.4.7) and the relation (3.4.19) is established. We remark also that the 
algebraic expression (3.4.11) for a2 is not necessary for the implementation of 
the algorithm. □ 



3.5 Entropy correction. 

• The Roe flux replaces the nonlinear waves of the gas dynamics, i.e. 
the rarefactions and the shock waves by linear waves that are the contact 
discontinuities. If sufficiently weak shock waves occur for a given discontinuity 
between two states VFieft and IKight, the Roe flux presented above is a good 
approximation, but if a rarefaction containing a sonic point is present among 
the nonlinear waves that solves the discontinuity problem between Wiqh and 
W^right, it has been early remarked that for this very particular situation, the 
Roe flux does not satisfy the entropy condition (see e.g. Godlewski and Raviart 
[GR96]). 

• A popular response has been proposed by Harten [Ha83] with a tuning 
parameter that plays in fact the role of an artificial viscosity and P. Roe himself 
[Roe85] has proposed a nonparameterized entropy correction for his flux. With 
G. Mehlman, we have treated the same subject by the introduction of hyperbolic 
nonlinear models with nonconvex flux functions and have proved a discrete en- 
tropy inequality if sufficiently weak nonlinear waves are present in the problem 
[DM96] . We detail here the modification of the algorithm that we have proposed 
and tested numerically for various gas dynamics problems. 
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• We introduce as above two states Wi = and Wr = and the 

Roe matrix A{Wi, Wr) described in the preceding sub-sections. We have in 
particular the relation 



(3.5.1) Wr-Wi = W^-W^ 



E 



:1 ^3 



and we do not make the confusion between the eigenvalues \j{W^) of the left 
state, XjiyV^) of the right state, and A* of the Roe matrix. We introduce the 

following two intermediate states and according to 

(3.5.2) = W^ + airl, = + ai rl -\- a2 r* = W^-a3r* 
and illustrated on Figure 3.2. Note that Aj(VF^) is well defined for j = 1, 2, 3 
and /c = 0, 1, 2, 3 : it is the eigenvalue of the intermediate state W'^. 
We define now the set S of sonic indices by the condition that the sign of the 

eigenvalue is increasing from negative to positive values accross some j-wave : 

(3.5.3) S = 2,,3}, Xj{W^-^) < < Xj{W^)}. 

The modification of the Roe flux is active only for the sonic indices and we 
introduce a polynomial pj of degree 3 by the classical Hermite interpolation 
conditions 

p,{o) = 0, p;.(o) = Xj{w^-') 

that defines explicitely the polynomial Pj{») by the algebraic relation 

XJW^) + XJW^-^)-2X* , 



(3.5.4) 



(3.5.5) 



+ 



3 A* 



2Xj{W^- 



■')-X,{W^) 



a. 



e + Xj{w^-') ^ 



K 


i 






\ 

\ 


/ y ^ 
1 y 

^ 3 









X 



Figure 3.2. 



Intermediate states for the entropy correction 
of the Roe upwind scheme. 



• With the hypothesis that j e S, it is not difficult to see [DM96] that 

the polynomial Pj{») has a unique minimum inside the interval (0, aj). The 
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argument of this point of minimum is given according to 
(3.5.6) Cj 



+ V (3 A* - X,{Wo) - X.iWo-^))' - X,{W^-^)X,{Wo) 
Since Pj{^*j) is the unique minimum of the polynomial Pj{») on the inter- 
val (0, aA, we have < and < A* . Then the modified flux 
$modif(^^^ VI^^) is defined from the Roe flux ^{Wu Wr) by the relation 

(3.5.7) $"^°^^^(iyi, Wr)^^{Wu Wr) + Y maxf ^illil, 

jes \ J 3 / 

that makes the added numerical viscosity explicit. 

3.6 Nonlinear flux boundary conditions. 

• At the two extremities a; = and a: = L of the pipe, we have to express 
on one hand the datum of a given nonstationary pressure n(t) at a: = and 
on the other hand a free output of the waves at x = L. 

• For the numerical boundary condition for pressure, we follow a general 
approach founded on the so-called partial Riemann problem [DuOl] that 
generalizes to nonlinear hyperbolic systems the reflection operator of relation 
(2.5.22). For a given discrete time = n At, and a given state = Wr in 
the first cell of the unidimensional mesh, we construct a boundary state Wq = 
Wi that satisfies the boundary constraint 

(3.6.1) p{W^) = n"+i/2, n>0, 

and moreover, we impose that the state W^^2 present in the first cell is issued 
from the boundary state Wq with an ingoing 3- wave, i.e. we impose the 
relation 

(3.6.2) W^/2 - = «3 

as illustrated on Figure 3.3. This problem has a unique solution, as claims the 

Proposition 3.7. Pressure flux boundary condition with Roe matrix. 

We consider a left boundary condition associated with a pressure 11 and a right 
datum defined by state Wr- Then there exists a unique left state Wi that 
satisfies the pressure condition 

(3.6.3) p{Wi) = n 

and such that when we construct the Roe intermediate state W* according to 
the relations (3.3.4) to (3.3.6), the difference Wr — Wi has only one component 
over the third eigenvector of the Roe matrix dF{W*), i.e. the relation (3.4.5) 
can been written under the form 
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^3 ^3 



(3.6.4) Wr - Wi 
The density pi and the velocity ui of the left state Wi are given according to 
the relations 

(7+i)n + {-f-l)Pr 



(3.6.5) 
(3.6.6) 



Pi 



Ui 



(7-i)n + (7+i)p, 

+ (n - Pr) 



Pr 



p,((7-l)n + (7 + l)Pr) 



U 



* * 

u - c 



Pressure 




* * 

u + c 



W 



1/2 



Figure 3.3. Nonlinear boundary condition for given pressure at x 

with the Roe upwind scheme. 







• The proof of Proposition 3.7 is a consequence precisely of the preced- 

ing subsections about the Roe flux. We first remark that the relations (3.4.5) and 
(3.6.4) are absolutly identical. Then we deduce that necessarily ai = a2 = ^ 
and according to the relations (3.4.11) and (3.4.12), we get 

(3.6.7) li — p* C* Ui = Pr — P* C* Ur 

(3.6.8) n-(c*)2pz = Pr-ic^^Pr 



We deduce simply Ur — ui 
(3.4.9), we get : 

(c*)2 = 



Pr 



n 



C* (Pr-Pl) 



p* C* 



and due to the relation 



P" 



1 



7-1 , f C*{pr-pi) \^ , f Pi + P* rj , P* + Pr 

-^p +7 — - — n+ Pr 



P' 



Pl 



pr 



Then after multiplication by {pr — Pi), we obtain with the help of (3.6.8) 



= (i?, - n) - ^ [Pr - n) ^ 



piY 



p* 

- VPI f Pl+P* TT , P* +Pr 

- 7 ^ , ^ n + Pr 



'Pr^ \fpl ^ 
7-1 

V - n) 



pl 

fPr 



+ 



'Pr 



pr 

- 7 



'Pr 



7 {Pr - n) 



-1 n - 7 1- 



'Pr 



Pr 
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We multiply the previous equality by . ^ and we get : 



2 pr p. 



id est 



7-1 



{Pr - n) + 7p, ^ = + ^ {pr - n) 



pr 

and the relation (3.6.5) is established 



2 

io: 

We deduce from the previous relation : 

(7+1) n + {^-i)vr 2(p,-n) 



^ - (7-i)n + (7+i)p. (7-i)n + (7+i)p. 

and due to the relation (3.6.8) : 

(^*^2 ^ (7-i)n + (7+i)p. ^ ^ / (7-i)n + (7+1) 



2pr 

and the relation (3.6.6) is an easy consequence of the last equality joined with 
(3.6.7). The proposition 3.7 is established. □ 

• The determination of a nonlinear nonreflecting boundary condition 

at x = L is still an open mathematical problem. We recommand for deriving a 
flux boundary condition for such a situation to impose that no wave are present 
at the interaction for the last interface j = J. We just write 

(3.6.9) /;+^/' = F(W-;_i/2) , n > 

which is equivalent of introducing a right boundary state Wj according to the 
simple relation Wj = Wj_-^^^2 then making these two states interacting 

with the Roe flux : /y+^Z^ = $(W'_y_i/2, W^). This last definition is equivalent 
to the one proposed in (3.6.9) due to the property (3.3.2) of the Roe matrix. 



4) Second order and two space dimensions. 
4.1 Towards second order accuracy. 

• The finite volume method described in the previous sections is a natural 
method for the discretization of systems of m conservation laws. It conducts 
to an explicit scheme in time : the evaluation of the field at time step 

(n+1) At needs only the knowledge of the field W^y^ for ji = 0, • • • , J— 1 
at the preceding time step nAt. This evaluation needs a certain number of 
auxiliary computations without the resolution of any linear system involving the 
new field. The method is parameterized by the choice of a numerical fiux and a 
great fiexibility can be adopted at this level. We have proposed two fiuxes for 
nonlinear problems related to nonlinear acoustics and gas dynamics, the Roe fiux 
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$(•, •) of relations (3.4.6)-(3.4.8) that conduct to a discrete scheme according 
to the relation 

(4.1.1) = $(w^;-i/2, , 

and the modified Roe fiux #) of relation (3.6.7) that enforces the 

entropy condition. This explicit version of the finite volume method is submitted 
to a stability condition that can be written as a first approximation for linear 
cases as : 

At 

(4.1.2) co^ < 1. 

• Nevertheless, the above finite volume method is only first order accu- 
rate. If we insert an exact solution W{x, t) of the conservation law 

dW d 

inside the formal expression of the fiux (4.1.1), it is easy to see that the finite 
difference ;^(/J+i^^^ — Z/^^^^) is first order accurate : 

(4-1.4) i^ifj:.''' - ir'^) = {c^)^Zl ^ ^ ■ 

In a similar way, the use of an explicit scheme in time conducts to 

(4.1.5) (w-j^^Y/ 2 - wjVv.) + ^ isTj' - fr") = 

and maintains this first order accuracy for the finite volume scheme. 

• We develop in this section the fact that it is possible to improve the 
method, i.e. to define a method with a relation of the type (4.1.5), and that 
conduct to a troncation error of second order : 

1 (yirn+l y^n ^ , 1 ( rn+1/2 ^n+l/2^ _ 



(4.1.6) { raw d \"+i/2 

The price to pay is to develop flux formulae much more complicated than the 
simple relation (4.1.1). When the second order precision (4.1.6) is achieved with 
a stable scheme, the precision is sufficient to develop predictive computations in 
acoustics and aerodynamics, whereas that is not the case with the initial scheme 
(4.1.1) (4.1.5). 



4.2 The method of lines. 

• The simplest way to extend the first order finite volume scheme is first to 
develop a new vision of the method with emphasis more on abstraction. We have 
presented a method founded on the integration of the conservation law (4.1.3) 
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inside the space-time domain V^^-^l^ = ]xj , Xjj^\[ x jt"^, t'^+^f as suggested in 

(2.1.6). With the method of hues, we just integrate the conservation (4.1.3) 
in space in each control volume = ]xj , Xjj^i[. It is straightforward to 

introduce the mean value Wj+i/2(t) in this finite element : 

1 r^j+i 

(4.2.1) Wj+y2{t) = y-^ 1/ W{x,t)&x- 

I I Jxj 

then we integrate the conservation law (4.1.3) in space in the cell Xj-|_i/2 and 
taking into account the relation ^Wjj^\/2{t) = Jx'^^ ^rC^j^)^^? we get 
simply 

(4.2.2) ^^.+i/2(i) + - Mt)] = 
with 

(4.2.3) f,{t) = F(W{xj,t)). 

• As usual with the finite volume method, a numerical scheme can be 
obtained from the relations (4.2.2) (4.2.3) by replacing the relation (4.2.3) by 
some explicit function over the set of all discrete variables introduced for the 
relation (4.2.1). To fix the ideas, we introduce a dynamic state vector Z{t) 
composed by all the dynamic variables on the finite mesh : 

(4.2.4) Z{t) = (T^i/2(t), •••,W^i+i/2(t), •••,H^j-i/2(t)) e OR^y. 
The discretization in space is achieved if we are able to determine the numerical 
fiux fj{t) with the help of both the dynamic state vector Z{») and the bound- 
ary conditions, that is the input pressure n(t) in the example considered in 
the last section. As in relation (2.1.12), we introduce a local numerical flux 
function ^j(», •) relative to the vertex xj : 

(4.2.5) f,{t) = ^,(U{t),Z{t)). 

We replace the relation (4.2.3) by the numerical approximation (4.2.5) inside 
the equation (4.2.2) of dynamic evolution of the state variable Wj+i/2(«). We 
obtain the following ordinary differential equation 

(4.2.6) ^w,+i/2{t) + ^[^,+i(n(t), z(t)) - *,(n(t),z(t))] = 0. 

• The method of lines is a semi-discrete version of the finite volume 
method. It is obtained by integration in space of the conservation law without 
integration in time. The result is not a numerical scheme but just an ordinary 
differential equation for the dynamic state vector Z{») described component by 
component with the equation (4.2.6). The method is parameterized by the local 
numerical flux functions •) and take the general form of a dynamical 
system parameterized by the pressure function 1 1 — > n(t) : 

(4.2.7) ^z(t) = G(Z, t) . 
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(4.2.12) 



The discrete dynamic function (IR^)-^ x [0, +00 [ 3 (Z, t) ^ G{Z, t) € 
^-^rny -g ^ vector valued expression with J components : 

(4.2.8) G{Z,t) = (Gi/2(Z,t), •••,G,-+i/2(Z,t), •••,Gj_i/2(Z,t)) e (IR^)^ 
and it is defined from the (J+1) local numerical fiuxes (^j)j=o, J with the 
very simple algebra relative to the finite volume method : 

(4.2.9) G,+i/2(^,t) = -^(^,+i(n(t), z)-^,(n(t),z)), o<j< 

Proposition 4.1. Explicit Euler scheme. 

With the choice of the first order scheme in space, that is 

(4.2.10) ^,(^^Z-) = $(T^;_i/2, W7+1/2) if j = !,•••, 

and the first order explicit forward Euler scheme for the ordinary differential 
equation (4.2.7), id est 

(4.2.11) (Z^+i - Z^) = G(Z", r) , 

we recover the previous first order finite volume scheme 

- *(W7-i/2> "^"+1/2)) =0 for i = 1, ■ ■ . , J-2. 

• We write the relation (4.2.10) for the particular control volume i^j+i/2 
and we get : ^7+i/2 = + ^iz"", due to (4.2.11), then 

^;+V2 = ^;+i/2 - - ^.■(^^ Z-)) due to (4.2.9) 

= w^;+i/2- ^(^(w^;Vi/2,w^A3/2)-^(^;-i/2,^;+i/^ c.f. (4.2.10) 

and the relation (4.2.12) is established. □ 
4.3 The method of Van Leer. 

• We turn now to the construction of a second order accurate version of the 
finite volume method as proposed initially with the "Multidimensional Upwind- 
centered Scheme for Conservation Laws" of B. Van Leer [VL79]. The fundamen- 
tal idea of this scheme is the reconstruction of a function IR 9 a: i — )■ W{x) € IR 
from his mean values in each cell • The reconstructed function 
is regular inside each control volume i^j+i/2 and is discontinuous at the in- 
terfaces xi between two control volumes. The application to the finite volume 
method replaces the scheme (4.1.1) by the same Roe fiux interaction $(•, •) 
considered for the two extrapolated data and Wj~ on each side of the 
boundary : 

(4.3.1) /, = ^wr,w;). 
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W(x,f) 




Figure 4.1. First order and second order interpolation at the interface xj. 

W(x,t") 





— • — 

w 

j-3/2 

®- 


T 
1 
1 

1 — • — 

1 w 

j+l/2 
<> 

j-1/2 ' 
1 


• 

W 

j+3/2 








1 

X. 




► 

X 



Figure 4.2. Construction of the nonhnear interpolated 
values WT^ and with the Van Leer method. 



• The simplest case is illustrated on Figure 4.1. It imposes simply the 

reconstructed function W{x) to be constant in each interval : 

(4.3.2) W{x) = 1^^+1/2 , Xj < X < Xj+i . 

The two limit values on each side of the interface xj are the following ones : 
= Wj_i/2 , Wj~ = and the explicit version of this finite volume 

scheme is the standard first order numerical fiux (4.1.1) as seen at Proposition 
4.1. In the context of the method of lines, we obtain : 

(4.3.3) = $(VF,-i/2, W^,+i/2). 
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• The second order accurate Muscl method consists first in restricting the 
methodology to a scalar field W{») and to construct an affine function in each 
interval Kj_^i^2 instead of a constant function as in (4.3.2). We set 

(4.3.4) W{x) = + Pj+i/2 {x - Xj^i/2) , Xj < X < Xj+i . 
The simplest choice for a slope is the one of the centered scheme : 

(4.3.5) p,+i/2 = ^(^,+3/2-1^,-1/2), 

and due to (4.3.4) and (4.3.5), the extrapolated values and Wj' on each 

side of the interface located at the position the following ones : 

(4.3.6) WJ- = Wj_^,2 + \ (T^, + l/2 - Vr,-3/2) 

(4.3.7) W+ = T^,-+i/2 - \ (T^,+3/2 - ^^,-1/2) . 

• The choice of a numerical flux given according to the relations 

(4.3.8) = ^(W,_y2 + \{W,+y2-W^-3/2).W,+y2-\{W^+3/2-^^ 

lead to an unstable scheme when we consider the particular case of the advec- 
tion equation with the first order explicit scheme in time. 

Proposition 4.2. Linear Muscl scheme is unstable. 

We apply the linear Muscl approach for the advection equation. Then the nu- 
merical scheme obtained by association of (4.3.8) and the upwind scheme (2.1.19) 
conducts to the following explicit first order scheme : 

This scheme is unstable for each At > 0. 

• Due to the expression (2.1.19) of the upwind scheme, the discrete first 
order in time advection equation can be written : 

w;vv2 - + ^ (w'^i" - w^-") = 

and the expression (4.3.9) is a consequence of the left extrapolation (4.3.6). For 
the study of stability, we introduce a profile of the type = ^ ^x) 

with a wave number k. The scheme (4.3.9) can be written as 

with an amplification coefficient ^(^, a) (^ = At, a = ^^) given simply by 
the expression 

g{^^(j) = 1 - a{l-e-'^) - ^(e^^-l-e-^^ + e-2^^), then 



(4.3.9) 
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^(^, cr) = 1 - (7(l-cosO - ^ (-l + cos2^) - i(7 (^sin^ + isin^-isin2^^ 
= l-(7(^l-cos^+i(2 cos^ ^ - 2 ) ) - i (7 sin^ - ^ sin^ cos ^ ^ , 



(4.3.10) g{^i,o) = 1 - -(1-cosO - i 2'i^^(3-^°'^)- 

For ^ arbitrarily small, we deduce from (4.3.10) : Idi^: <^) 1^ = 1 + cr^ + 
0(^^) which establishes the instability for all cr ^ 0. □ 



• The above remark motivates the introduction of so-called slope lim- 

iters, intensively studied during the period 1980-90 after the pioneering work 
of Van Leer [VL77]. The idea is to search an interpolation of the field 

W{») at the left of the point xj from the neighbouring mean values Wj-3/2, 
Wj_i/2 and and by left-right invariance of the procedure, to con- 



struct an interpolated value Wj from the first right neighbours Wj_i/2, 
Wj+1/2 and Wj+3/2 as suggested on Figure 4.2. We replace the relations 
(4.3.6) and (4.3.7) by a nonlinear interpolation parameterized by a slope 
limiter IR 9 r 1 — > (p{r) G IR : 



(4.3.11) W. 



1 



(4.3.12) 



J 



-1/2 + 2 ^' 
= Wj+i/2 - 



w. 



3+^/2 



w. 



J -1/2. 







'3 + 1/2 - 

We remark that the relations (4.3.6) and (4.3.7) are a particular case of the gen- 
eral nonlinear relations (4.3.11) and (4.3.12) with the particular choice ^{r) = 
I (1 -|- r) . The limiter function satisfies very often the functional relation 



- W-,-1/2) 



(4.3.13) ^{r) = r(f(^^^ , r > 0. 



Among all the possible choices, we have adopted for fluid mechanics [DM92] the 
so-called STS-limiter deflned by the relations 

, 
3 

T' 

1 + r 



(4.3.14) 



r < 
< r 



3 
2 



1 

r > 2 



1 

< - 
- 2 



< 2 



and illustrated on Figure 4.3. 
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Figure 4.3. Examples of limiter functions that can be easily extended 

to unstructured meshes. 



4.4 Second order accurate finite volume method for fiuid problems. 

• We detail in this section a generalization for unstructured meshes of 

the Muscl scheme proposed by Van Leer [VL79]. At one space dimension on a 
uniform mesh, it is classical to consider a scalar field z among the primitive 
variables, id est 

(4.4.1) z e (primitive variables) 

and instead of computing the interface flux with relation (4.1.1), to first construct 
two interface states Wg and Wg on each side of the interface S. Then the 
flux is evaluated by the decomposition of the discontinuity : 

(4.4.2) fs = ^{Ws ,W+), Se{xu---, xj_i } . 

This nonlinear interpolation is done with a slope limiter (/?(•) that operates 
on each variable proposed in (4.4.1) and we have typically when a left-right 
invar iance is assumed [Du91] : 

(4.4.3) zg = ^,-1/2 + fe-+i/2-^j--i/2) , S 

(4.4.4) 4 = ^J^U2-y Cr''~T''' ){^,^U2-^J-i/^), S- 



Xj 



Xj 
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Figure 4.4. Structured Cartesian mesh. 
The control volumes are exactly the elements of mesh T. 




Figure 4.5. Unstructured mesh composed by triangular elements. 
The control volumes are exactly the elements of the mesh. 



• We focus now on the use of unstructured meshes for the extension to 

second order accuracy of the finite volume method. As in the one-dimensional 
case, the domain of study is decomposed into finite elements (or control volumes) 
K G Sj- than can be structured in a Cartesian way (Figure 4.4) or with a cellular 
complex as in Figure 4.5. In both cases, the intersection of two finite elements 
define an interface f ^Tj-- We denote by n/ the normal at the interface / that 
separates a left control volume Ki{f) and a right control volume Kr{f). The 
ordinary differential equation (4.2.6) is replaced by a multidimensional version : 

(4.4.5) \K\ ^+ Yl \f\ ^(WK,nf,WK^^f^) = 0, K e Sr ■ 

fCdK 
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For internal interfaces, the function $(•, nj , •) is equal e.g. to the Roe flux 
between states and in the one-dimensional direction along nor- 

mal in order to take into account the invariance by rotation of the equations 
of gas dynamics (see [GR96]). 

• We consider now a finite element K internal to the domain. The exten- 
sion to second order accuracy of the finite volume scheme consists in replacing 
the arguments and in relation (4.4.5) by nonlinear extrapola- 

tions Wxi{f),f and Wx^{f),f on each side of the boundary of state data and 
evaluated as described in what follows. We first introduce the set J\f{K) of 
neighbouring cells of given finite element K e £-r, as illustrated on Figure 4.6 : 

(4.4.6) Af{K) = {LeSr, eTr, f c dKndL} . 

For L e J\f{K), we suppose by convention that the normal Uf to the face 
/ C OK n dL is external to the element K id est Kr{f) = K, Ki{f) = L. 
We introduce also the point yK,f on the interface / C dK that links the 
barycenters xk and Xx^{f) '■ 

(4.4.7) I yK,f = {'^-OK,f)xK + OKjXK^^f), VKJ^f, 

y / C dK , K finite element internal to mesh T . 

Then, following Pollet [Po88], for z equal to one scalar variable of the family : 

(4.4.8) z e {p, pu, pv , p} 

we evaluate a mean value zk, / on the interface / : 

(4.4.9) ZWJ = {l-eK,f)zK + eK,fZK^^f) 

and the gradient \/z{K) of field z{») in volume K with a Green formula : 




Figure 4.6. Cellular complex mesh with triangles and quadrangles. 
Three neighbouring cells are necessary to determine the gradient 
in triangle K and to limit eventually its variation. 
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• An ideal extrapolation of field z{») at the interface / would be : 

(4.4.11) ZK,f = zk + Vz{K).(yK,f-XK) 

but the corresponding scheme is unstable as explicited at Proposition 4.2. When 
the variation 'S/z{K)» {yK,f — xk) is very important, it has to be "limited" 
as first suggested by Van Leer [VL77] . For doing this in a very general way, we 
introduce the minimum raK{z) and the maximum Mk{z) of field z{») in the 
neighbouring cells : 

(4.4.12) rriKiz) = min { , LeX{K)} 

(4.4.13) Mk{z) = max{zL, LeAf{K)}. 

If the value zk is extremum among the neighbouring ones, ie. if zk < rriKiz) , 
or Zk > Mk{z), we impose that the interpolated value ZK,f is equal to the 
cell value zk '■ 

(4.4.14) ZK,f = Zk if zk < mK{z) or zk > Mk{z) , f C dK . 
When on the contrary zk lies inside the interval [mK{z), Mk{z)], we impose 
that the variation zk, f — zk is limited by some coefficient k {0 < k < 1) 
multiplied by the variations zk — rnK{z) and Mk{z) — zk- We introduce a 
nonlinear extrapolation of the field ^(•) between center xk and boundary face 
VK, f if C dK) : 

(4.4.15) ZK,f = ZK + aK{z)Vz{K).(yK,f-XK), f ^ dK 
with a limiting coefficient aK{z) satisfying the following conditions : 

{0 < aK{z) < 1 , z{*) scalar field defined in (4.4.8), K e Er 
k{zK-mK{z)) < aK{z)Vz{K).(yK,f-XK) < k{MK{z)-ZK) 
\ffcdK, KeSr- 
Then aK{z) is chosen as large as possible and less than or equal to 1 in order 
to satisfy the constraints (4.4.16) : 

m.m{MK{z) - Zk, Zk - mK{z) ) 



(4.4.17) aK{z) = min 



1, k 



max { I Vz{K) • [yK, f - xk) \ , f ^ dK] 



• In the one dimensional case with a regular mesh, it is an exercice to 

re- write the extrapolation (4.4.15) under the usual form (4.4.3) in the context 
of finite differences. In this particular case, some limiter functions r i — )■ (pk{r) 
associated with particular parameters k are shown on Figure 4.3. For k = 1, 
we recover the initial limiter proposed by Van Leer in the fourth paper of the 
family "Towards the ultimate finite difference scheme..." [VL77] ; for this reason, 
we have named it the "Towards 4" limiter (see Figure 4.3). When k = ^ we 
obtain the "min-mod" limiter proposed by Harten [Ha83]. The intermediate 
value k = ^ is a good compromise between the "nearly unstable" choice k = 1 
and the "too compressive" min-mod choice. We have named it STS (see also 
(4.3.14)) and it has been chosen for our Euler computations in [DM92]. 
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Figure 4.7. Slope limitation at a fluid boundary. 




Figure 4.8. Slope limitation at a solid boundary. 

• We explain now the way the preceding scheme is adapted near the bound- 
ary. We first consider a fluid boundary. When K is a, finite element with some 
face g C dK lying on the boundary, we still define the set J\f{K) of neigh- 
bouring cells by the relation (4.4.6) as shown on Figure 4.7. The number of 
neighbouring cells is just less important in this case. Then points Vk,/ are 
introduced by relation (4.4.7) if face / does not lie on the boundary and by 
taking the barycenter of face g if it is lying on the boundary. The only differ- 
ence is the way the values zk, g are extrapolated for the face g that is on the 
boundary ; we set 

= Zk , g cdK , 

g face lying on the boundary of the domain. 
When values zk, / are determined for all the faces / C dK, the gradient 




(4.4.18) 
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Vz{K)^ the minimal mK{z) and maximal Mk{z) values among the neigh- 
bouring cells are still determined with the relations (4.4.10), (4.4.12) and (4.4.13) 
respectively. The constraints (4.4.16) remain unchanged except that no limita- 
tion process is due to the faces lying on the boundary. In a precise way, we 
set : 

{axiz) = 
^.^L k mm{MK{z) - Zk , Zk - mK{z)) ' 
""'""L ' ma^ {\Vz{K).{yK, f-XK)\,fC OK, Kr{f) G Af{K)}\ ' 
Then the interpolated values zk, / for all the faces / C dK are again predicted 
with the help of relation (4.4.15). 

• For a rigid wall, the limitation process is a little modified, as presented 
at Figure 4.8. We first introduce the limit face g inside the set of neighbouring 
cells ' 

(4.4.20) A.W= {^^^-^/^f . H ; 

[g^J-j-, g C oK, g on the boundaryj. 

For the face(s) g C dK lying on the solid boundary, we determine preliminary 
values Zk, g by taking in consideration at this level the impenetrability bound- 
ary condition u • — 0. We introduce the two components Ug and of 
the normal at the boundary and we set, in coherence with variables (4.4.8) : 

i PK^ = PK 
PK,g UK,g = PK {uk - {Uk •Tig) u'^g) 
PK^ vic^ = PK (vK - {uk • rig) ) 
PK,g = Pk- 

We consider also these values for the limitation algorithm. We define "external 
values" zl for L = g and face g lying on the boundary as equal to the ones 
defined in relation (4.4.21) : 

(4.4.22) Zg = Zk, g , z{») field defined in (4.4.21), c c^iT on the boundary. 
Then the extrapolation algorithm that conducts to relation (4.4.15) for extrap- 
olated values Zk, / is used as in the internal case. 

4.5 Explicit Runge-Kutta integration with respect to time. 

• When all values zk, / are known for all control volumes K e ^7-, all 
faces f dK and all fields z{») defined at relation (4.4.8), extrapolated states 
Wk, f are naturally defined by going back to the conservative variables. Then 
we introduce these states as arguments of the flux function $(• , ny , •) and 
obtain by this way a new system of ordinary differential equations : 

(4.5.1) \K\^+ Yl \f\^(WK,f:nf,WK^^f),f) = 0, Keer- 

fCdK 
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• The numerical integration of such kind of system is done with a Runge- 

Kutta scheme as presented in [CDV92] . We have used with success in [DM92] the 
Heun scheme of second order accuracy for discrete integration of (4.5.1) between 
time steps nAt and (n+l)At : 

(4.5.2) ^(^-W^)+E Ifl^i^lf^^f^^KM),/) =^^^^^r 

fCdK 

(4.5.3) \^ I^Wk-Wk)+ Yl \f\^{wK,f:nf,WK^^f^,f) =0,Keer 

fCdK 

(4.5.4) W^+' = 1(^Wk + W^'), KeSr- 
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